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1.  Objective 


Better  understanding  the  stable  boundary  layer  (SBL)  through  high-resolution  numerical 
modeling  is  the  goal  of  this  report.  At  all  phases,  the  model  results  are  compared  to  analytical 
results  and  actual  observations.  The  numerical  computations  should  reveal  SBL  details, 
requiring  further  experimental  verification  and  model  enhancements. 


2.  Approach 


Currently,  SBL  processes  are  simulated  to  provide  a  virtual  laboratory  to  develop  quantitative 
predictions  on  the  role  of  various  processes  commonly  observed  in  stable  environments. 

Specific  focus  is  on  how  gravity  waves,  generated  by  shear  and  terrain,  may  alter  the  state  of 
localized  regions  within  the  SBL,  thus  allowing  turbulence  to  develop.  Simulation  results  are 
being  carefully  compared  with  current  observational  data,  other  numerical  treatments,  and 
laboratory  experiments  to  refine  modeled  terms.  Successful  simulation  should  stimulate  new 
experimental  efforts  in  verification  and  model  refinement.  Contrary  to  general  perception, 
temperature  and  wind  profiles  drawn  through  relatively  coarse  resolution  observations  may  belie 
the  presence  of  small-scale  structures  as  seen  in  the  acoustic  sounder  observations  in  figure  1 . 


Figure  1.  Acoustic  sounder  profiles  showing  turbulent  activity  deformed  gravity  waves  (7). 

Observations  have  shown  sublayers  where  temperature  and  wind  gradients  are  more  favorable 
for  instability  (figure  2).  These  layers  are  then  deformed  by  mesoscale  motions  propagating 
through  the  SBL  (2).  How  these  features  fonn  is  poorly  understood,  but  gravity  waves 
(ubiquitous  in  the  SBL)  may  be  one  responsible  mechanism.  During  the  month-long  CASES-99 
experiment,  wave  activity  was  observed  nightly  (J).  Gravity  waves  transporting  momentum  and 
energy  without  loss  is  well  documented  (2,  4 ,  5).  What  is  investigated  here  is  whether  gravity 
waves,  created  internally  or  externally  to  the  SBL,  can  alter  the  state  of  small  regions  within  the 
boundary  layer,  fostering  instability.  The  resulting  turbulence  becomes  part  of  a  complicated 
system  of  energy  and  momentum  exchange  between  the  turbulence,  waves,  and  mean  flow. 


1 


Source:  Boulder  Atmospheric  Observatory. 

Figure  2.  Observations  using  a  translatable 

boom  with  meteorological  instruments 
(2). 

If  the  energy  and  momentum  transfer  is  between  the  SBL  and  free  atmosphere,  universal  scaling 
will  not  be  very  accurate  (6). 

Current  efforts  center  on  describing  the  fine-structure  evolution  instead  of  improving  bulk  SBL 
parameterizations.  Initial  focus  is  on  the  shear  instability  and  terrain  effects  that  produce  gravity 
waves.  To  explore  wave-turbulence  interactions,  the  model  must  resolve  small-scale  structures 
yet  properly  generate  and  trap  gravity  waves.  Mechanisms  for  the  ducting  of  gravity  waves 
include  critical  levels  (5),  Doppler  ducting  (7),  and  thermal  ducting  (if  the  SBL  is  bounded  by 
neutral  or  convective  lapse  rates).  The  National  Taiwan  University/Purdue  University 
(NTU/PU)  nonhydrostatic  model  (8,  9),  modified  to  include  a  more  complete  treatment  of 
turbulence  in  stably  stratified  environments,  was  used.  The  model  has  been  developed  for  more 
than  10  years  and  applied  to  several  different  physical  situations,  such  as  nonlinear  mountain 
waves,  shallow  convection  due  to  cold-air  outbreaks,  deep-vortex  generation,  and  shedding. 

The  NTU/PU  model  uses  the  prognostic  equation  for  density  instead  of  a  pressure  tendency 
equation.  The  method  also  explicitly  resolves  acoustic  waves  without  time  filtering.  Due  to  the 
small  time  step  required  to  satisfy  the  Courant-Friedrichs-Lewy  criterion  for  acoustic  waves,  two 
methods  are  used  to  conserve  computational  resources.  By  decreasing  the  speed  of  sound,  the 
largest  stable  time  interval  can  be  increased  without  affecting  the  meteorological  evolution. 

Also,  the  model  uses  a  split  time  step  integration  scheme  to  update  slower  varying  physics  in  a 
more  appropriate  interval.  Finally,  by  avoiding  the  difficulties  of  solving  pressure  tendency 
equations,  the  model  is  more  easily  parallelized  for  efficient  use  of  massively  parallel  clusters  (5) 
and  uses  a  stretchable  (in  all  three  directions)  terrain-following  grid.  Using  a  nested  high- 
resolution  region  allows  more  complicated  lateral  and  upper  boundary  conditions  to  be  specified, 
which  is  especially  important  when  simulating  conditions  to  prevent  gravity  waves  from 
propagating  out  of  the  system.  The  terrain-following  vertical  coordinate  (g-z)  is  prescribed  at 
model  initialization  to  include  terrain  effects. 


2 


In  environments  with  strong  stratification,  turbulent  motions  are  damped  in  the  vertical  direction 
leading  to  anisotropy  of  the  turbulent  field  (10).  In  an  eddy-viscosity  closure  scheme,  flow 
history  takes  a  reduced  role;  fluctuations  are  mostly  dependent  on  the  instantaneous  mean  strain 
rates.  Isotropic  flow  history  only  enters  through  the  variables  used  to  calculate  eddy  viscosity 
and  eddy  diffusivity,  usually  by  turbulent  kinetic  energy  (TKE)  and  dissipation.  NTU/PU’s 
default  turbulence  scheme  model  is  a  single  equation  eddy-viscosity  scheme.  TKE  is  solved 
prognostically  while  the  length  scale  is  prescribed  and  is  based  on  model  grid  spacing  and 
stratification.  It  is  now  computationally  feasible  to  prognostically  calculate  transport  equations 
for  turbulent  moments  (a  method  developed  in  the  70s  by  several  groups).  The  method 
employed  in  the  NTU/PU  model  was  developed  by  Lumley  (11-13).  Several  terms  of  the 
second-order  equations  can  be  calculated  explicitly.  Some  of  the  other  terms — deviatoric 
dissipation,  isotropic  pressure  dissipation,  and  molecular  diffusion — are  neglected.  Deviatoric 
pressure  dissipation  is  modeled  using  a  retum-to-isotropy  and  rapid  strain  with  buoyancy 
contributions.  The  second-order  model  equations  are  as  follows: 


the  TKE: 

S,e  +  Ufdte  =  -u'lu\Sit  +  u\ff  -  e  ; 

the  potential  temperature  variance: 

a,¥+u,Sl¥-^- 


(1) 


(2) 


the  heat  flux: 


dt  u’jO’  +  U(8t  u'jff  =  -—u'(ffSej  - 


f  2  1 

btj  +  -eSf. 
v  J  3  V 


J  r.  7  3  j30n 
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and  the  deviatoric  momentum  flux: 


d,bu  +  U.dA.  =  ~—eS„  -  -  Z,  -  R„  -  >  -  ^-b„  +  —B„  ; 

1  v  z  z  tj  \  5  lJ  5  lJ  lJ  lJ  t  ^  10  J 
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(5) 


This  method  is  augmented  by  the  dissipation  spectra  for  stable  stratification  developed  by 
Canuto  and  Minotti  (10).  For  coarse  resolution  simulations  under  strong  stratification,  a  volume 
element  is  too  large  to  resolve  isotropic  turbulent  fluctuations.  Applying  the  Kolmogorov 


3 


similarity  for  the  inertial  subrange  overestimates  the  amount  of  energy  cascading  to  smaller 
scales;  the  subgrid  turbulent  motions  work  against  the  stratification,  leaving  less  to  cascade.  As 
grid  spacing  becomes  much  smaller  than  the  buoyancy  length  scale  ( Ab  =  Ve/A),  unresolved 
motions  become  isotropic  and  dissipation  again  follows  the  Kolmogorov  spectrum.  The 
functional  form  of  the  dissipation  is  as  follows: 

e3/2 

s=C,(e)—  ,  (6) 

and 

— -e1/2 

^C^iey2  —  ,  (7) 


with  dissipation  coefficients, 


Cs(e)  =  cEexA 


-0.053 


a2a2 


and  Ce6  (e)  =  C; 
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V  01. 


L Ce 
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a/a2 
e  ) 
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,  and  A  =min(Ax,Av,Az). 


(9) 


The  derivation  of  the  second-order  transport  equations  does  not  introduce  any  new  information 
into  the  equation  system;  previously  implicit  forcing  is  simply  exposed.  As  such,  closure  is  not 
introduced  due  to  the  presence  of  the  third-order  moments,  the  turbulent  transport  terms.  A 
similar  process  yields  third-order  transport  equations.  To  introduce  closure,  a  quasi-Gaussian 
assumption  is  employed  to  rewrite  the  fourth-order  moments  as  a  function  of  the  second-order 
moments  (11).  Over  the  bulk  of  the  SBL,  these  terms  should  have  a  small  effect;  however,  in 
regions  of  instability,  the  third-order  tenns  may  possibly  influence  the  evolution  of  the  small- 
scale  structures.  Due  to  the  computational  requirements  for  prognostically  solving  such 
equations,  preliminary  simulations  have  completely  neglected  third-order  correlations.  However, 
given  the  recently  completed  parallel  version  of  the  turbulence  closure  module,  third-order 
equations  can  now  be  practically  integrated. 

Although  initial  feasibility  testing  of  the  three-dimensional  higher-order  closure  model  (14)  was 
done  in  the  serial  version  of  the  NTU/PU  model,  sensitivity  testing  and  probing  was  not  practical 
given  the  serial  model’s  limited  resources.  The  new  closure  module  has  been  rewritten  to 
employ  the  message-passing  interface  (MPI)  libraries  for  parallel  execution  and  has  been  added 
to  the  parallel  version  of  the  NTU/PU  model.  The  parallel  version  enables  the  large 
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computational  domains  and  high  resolutions  required  to  capture  the  larger  and  smaller  scale 
processes  observed  in  the  SBL.  So  far,  the  performance  numbers  are  favorable  and  should  not 
hinder  future  large  simulations.  An  example  is  a  domain  with  1200  x  150  horizontal  grid  points 
and  60  vertical  levels.  Figure  3a  shows  the  performance  gained  by  increasing  the  number  of 
processors  leaving  the  domain  size  constant. 


The  better-than-double  performance  is  likely  attributed  to  better  processor  cache  usage  with 
smaller  domains  and  smaller  amounts  of  data  copied  in  the  custom  variable  packing  routines.  Of 
course,  for  increasing  number  of  processors,  CPU  resources  are  offset  by  an  increase  in 
communication  overhead,  thus  reducing  performance.  Since  our  focus  is  on  running  high- 
resolution,  large-domain  simulations,  a  more  relevant  comparison  is  to  test  scalability  by 
choosing  a  reasonable  subdomain  for  each  processor.  Therefore,  as  a  greater  number  of 
processors  are  used,  the  domain  increases.  This  prevents  the  total  domain  from  being  divided 
into  pieces  too  small  to  be  efficient  as  the  number  of  processors  becomes  large.  Hence,  the 
domain  on  each  processor  was  held  at  30  x  25  horizontal  grid  points  with  60  vertical  levels. 
Ideally,  as  the  number  of  processors  increases,  the  computational  time  should  remain  constant; 
however,  the  overhead  needed  for  parallel  communication  prevents  the  ideal  case.  As  figure  3b 
shows,  employing  960  processors  leads  to  an  1 1%  decrease  in  efficiency,  neglecting  output.  The 
final  simulation  domain  was  1200  x  600  x  60  grid  points.  Data  input/output  (I/O)  also  required 
careful  consideration.  I/O  was  treated  using  the  Hierarchical  Data  Format  v5  (HDF5)  libraries 
developed  at  the  University  of  Illinois  National  Center  for  Supercomputing  Applications.  HDF5 
is  an  abstracted  I/O  library  that  can  create  structured,  self-describing,  random-access  files  using 
independent  and  collective  (parallel)  I/O  operations. 
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Figure  3.  Input  performance  results  for  (a)  constant  total  domain  split  over  an  increasing  number  of 
processors  and  (b)  a  domain  that  increases  proportional  to  the  number  of  processors. 


3.  Results 


The  initial  simulations  used  a  serial  model  version  (14);  extensive  testing  of  the  closure  module 
was  impractical.  More  testing  is  needed  to  confirm  that  the  results  are  physical  and  not  caused 
by  numerical  issues  (e.g.,  contamination  by  lateral  boundary  conditions).  Kelvin-Helmhotz 
(KH)  instability  is  the  dominant  SBL  turbulence  mechanism;  it  results  from  wind  shear  present 
from  the  bottom  of  the  free  atmosphere  extending  into  the  boundary.  Thin  layers  of  relatively 
strong  shear  may  develop  as  gravity  waves  dump  momentum  into  localized  regions;  thus, 
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properly  simulating  shear  instability  is  a  primary  concern.  The  simulations  were  initialized  with 
a  100-m-deep  shear  layer  in  the  middle  of  the  domain,  with  a  4.9  m/s  velocity  difference.  The 
shear  region  velocity  profile  was  defined  by  a  hyperbolic  tangent  function;  the  thermal  gradient 
was  defined  by  a  constant  Brunt-Vaisala  frequency.  A  random  perturbation  was  added  to  the 
initial  u-wind  field.  Using  the  single-equation  eddy-viscosity  closure  scheme,  the  turbulent 
diffusion  smoothes  the  perturbations,  preventing  amplification  of  the  preferred  mode  and 
resulting  in  TKE  profiles  (figure  4a  and  b). 


Figure  4.  TKE  profiles  simulated  using  an  eddy-viscosity  closure  for  (a)  35 
and  (b)  48  min.  Potential  temperature  (to  show  overturning)  using 
second-order  closure  for  (c)  25  and  (d)  35  min. 

Switching  to  a  second-order  closure  reveals  a  different  evolution  (figure  4c  and  d).  By  25  min 
into  the  simulation,  a  narrow  band  of  modes  is  clearly  present.  By  35  min,  the  waves  are 
breaking,  forming  KH  billows.  A  spectral  decomposition  shows  the  energy  distribution  during 
the  wave  growth  and  breaking  phases  (figure  5).  Linear  theory  predicts  the  amplification  of  a 
single  preferred  mode  proportional  to  the  shear  layer  depth  (the  red  line  on  the  spectral  plots) 
(15).  In  the  initial  wave  growth,  a  narrow  band  of  modes  surrounding  the  linear  mode  is 
stimulated.  Later,  during  wave  breaking,  the  linear  mode  is  still  dominant,  but  secondary  modes 
also  appear,  especially  at  lower  wave  numbers,  which  have  been  observed  in  nature  (figure  6). 
Breaking  waves  are  the  primary  instability  associated  with  shear;  however,  several  other  effects 
have  been  observed:  subharmonic  vortex  pairing,  convective  rolls  development  with  vorticity 
perpendicular  to  the  vorticity  of  the  billow,  vortex  knotting,  and  vortex  tube  formation  (16). 
With  the  new  MPI  model,  simulating  such  effects  to  support  the  model’s  ability  to  accurately 
treat  subtle  processes  will  be  attempted. 
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Figure  5.  Fourier  modes:  (a)  25  min  during  initial  wave 
excitation  and  (b)  35  min  during  wave  breaking. 


Figure  6.  Frequency  modulated  continuous  wave  radar  image  of  KH 
waves  showing  the  presence  of  two  separate  modes  (2). 


4.  Conclusions 


The  SBL  continues  to  elude  quantitative  description;  much  work  remains  for  theorists, 
experimentalists,  and  modelers.  Given  the  observed  fine-scale  structures,  numerical  modeling 
may  provide  a  cost-effective  way  of  probing  the  behavior.  The  overall  goal  of  this  project  is  to 
simulate  this  complicated  behavior.  The  interactions  between  line-scale  turbulence,  larger-scale 
features  provided  the  impetus  for  implementing  a  new  higher-order  turbulence  closure  routine 
into  the  NTU/PU  model.  Preliminary  results  from  the  series  version  of  the  code  were 
encouraging,  but  the  inability  to  take  advantage  of  massively  parallel  clusters  significantly 
limited  the  scope  of  testing.  An  efficient  parallel  version  of  the  model  code  has  been  developed. 
Using  the  U.S.  Army  Research  Laboratory  HPC  clusters,  expanded  simulations  of  shear 
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instability  supported  the  utility  of  the  new  scheme.  The  model  showed  rapid  amplification  of 
perturbations  in  a  narrow  band  of  preferred  nodes  in  agreement  with  linear  theory.  The 
amplification  led  to  the  development  of  wave  breaking,  Kelvin-Helmholtz  billows  and  the 
eventual  decay  of  the  turbulence.  Because  shear  instability  is  the  primary  source  of  turbulence  in 
the  SBL  and  the  most  common  mechanism  for  the  production  of  gravity  waves,  accurate 
simulation  of  shear  instability  will  play  an  important  role  in  probing  gravity  wave-turbulence 
interactions. 

Much  work  remains  to  be  done,  and  more  comprehensive  simulations  are  being  planned.  An 
important  goal  is  to  demonstrate  the  destabilization  of  thin  sublayers  embedded  within  the  SBL 
via  in-bound  gravity  waves.  The  effects  of  nonuniform  terrain  and  terrain-generated  gravity 
waves  will  also  be  included.  Later  in  2008,  more  holistic  runs  will  be  attempted  by  initializing 
the  model  with  observations  from  CASES-99  field  experiments.  These  model  runs  will  provide 
an  important  test  of  the  utility  of  the  model  in  more  realistic  environments  and  will  hopefully 
produce  quantitative,  testable  predictions  that  will  facilitate  new  theoretical  work  on  transient 
turbulence,  turbulence  under  stable  conditions,  and  new  field  experiments.  Ultimately,  we  hope 
to  contribute  to  better  forecasts  of  meteorological  conditions  at  night,  when  the  Army  prefers  to 
operate,  and  the  dispersion  and  diffusion  of  chemical,  biological,  or  radiological  agents. 
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